sampleSize=1000;
outputSize=1000000;
lagMax=4; #Max lag for autocor
ro<-seq(length=lagMax+1, 0, 0); #autocor
Z<-seq(length=outputSize, 0, 0); #output standard gaussian
Y<-seq(length=lagMax, 0, 0);

#temp
Y[1]=1;
Y[2]=1;
Y[3]=1;
Y[4]=1;

#build covariance matrix (inversely), M means N-1
ro[1]=1;
ro[2]=0.6;
ro[3]=0.8;
ro[4]=0.4;
ro[5]=0.5;

#build covariance matrix (inversely)
Ro_N_N=seq(length=(lagMax+1)*(lagMax+1), 0, 0);
attr(Ro_N_N, "dim")=c(lagMax+1, lagMax+1);
for(i in seq(1,lagMax+1,1))
{
	for(j in seq(1, lagMax+1, 1))
	{
		Ro_N_N[i, j]=ro[abs(j-i)+1];
	}
}
Ro_1_2=Ro_N_N[c(1:lagMax),c((lagMax+1):(lagMax+1))];
attr(Ro_1_2, "dim")=c(lagMax,1);
Ro_2_1=Ro_N_N[c((lagMax+1):(lagMax+1)), c(1:lagMax)];
attr(Ro_2_1, "dim")=c(1,lagMax);
Ro_2_2=Ro_N_N[c((lagMax+1):(lagMax+1)), c((lagMax+1):(lagMax+1))];
Ro_1_1=Ro_N_N[c(1:lagMax), c(1:lagMax)];

#initialize output
mu_matrix=seq(length=(lagMax),0,0);
for(i in seq(1, lagMax, 1))
{
	Z[i]=Y[i];
	mu_matrix[i]=Y[i];
}
attr(mu_matrix, "dim")=c(lagMax,1);

temp1=Ro_2_1 %*% (solve(Ro_1_1));
new_var=Ro_2_2-(Ro_2_1 %*% (solve(Ro_1_1)) %*% Ro_1_2);
var_error=0;
if(new_var<=0)
{
	var_error=1;
}
new_ro=sqrt(new_var);

#calculate new mu and rio
for(i in seq(lagMax+1, outputSize, 1))
{
	new_mu=temp1 %*% mu_matrix;
	Z[i]=rnorm(1, mean=new_mu, sd=new_ro);
	
	for(j in seq(1, lagMax, 1))
	{
		mu_matrix[j,1]=Z[i-lagMax+j];
	}
}

#calculate autocor
Z_acf=acf(Z, lag.max=lagMax, type="correlation", plot=FALSE);
print(Z_acf);

#trytry based some func_y_to_x
W=func_y_to_x(Z);
trace1=W;